﻿using LxBasic;
using System;

namespace LxCore
{
    public class EquationWithDiffTerm : EquationWithTransTerm
    {
        public IField<double> Gamma;
        public void SetZone_DiffTerm(IField<double> diffTerm)
        {
            this.Gamma = diffTerm;
        }
        protected void formABOnFixValueWallBound(FvBound bound, FixValueWallBound fixValueWallBound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var gammaf = Gamma.GetFaceVarArr(bound);
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                double gamma = gammaf.GetFaceVar(i1);
                double xbou = BouX[i1];
                double alun = ABouH[i1] = -gamma * face.gDiff;
                BouFaceShare[i1] = gamma * face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, xbou) - alun * xbou;
            }
        }
        protected void formABOnFluxWallBound(FvBound bound, FluxWallBound fluxBound, ThreadInfo threadInfo)
        {
            //formABOnFluxWallBound_Test(bound, fluxBound, threadInfo); return;
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var gf = Gamma.GetFaceVarArr(bound);
            int start = bound.FaceIdx.Start;
            if (fluxBound.IsConstSource())
            {
                var factors = fluxBound.BouFactors1;
                if (fluxBound.FixFlux_b.IsEqual0())
                {
                    Array.Clear(ABouH, i1, endp1 - i1);
                    Array.Clear(BouFaceShare, i1, endp1 - i1);
                    for (; i1 < endp1; i1++)
                    {
                        ref var face = ref BouFaces[i1];
                        double beforPhib = face.gDiff;
                        double rightTerm = -face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, BouX[i1]);
                        double B = rightTerm / beforPhib;
                        factors[i1 - start] = (B);
                    }
                }
                else
                {
                    Array.Clear(ABouH, i1, endp1 - i1);
                    var fb = fluxBound.FixFlux_b.GetFaceVarArr(bound);
                    for (; i1 < endp1; i1++)
                    {
                        ref var face = ref BouFaces[i1];
                        double gamma = gf.GetFaceVar(i1);
                        double flux_b = fb.GetFaceVar(i1);
                        flux_b *= face.Area;
                        double beforPhib = face.gDiff * gamma;
                        double rightTerm = -gamma * face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, BouX[i1]) + flux_b;
                        double B = rightTerm / beforPhib;
                        factors[i1 - start] = (B);//phi_b=A*phi_C+B
                        BouFaceShare[i1] = flux_b;
                    }
                }
            }
            else
            {
                var factors = fluxBound.BouFactors2;
                var fb = fluxBound.FixFlux_b.GetFaceVarArr(bound);
                var fa = fluxBound.FixFlux_a.GetFaceVarArr(bound);
                var fh = fluxBound.Fixh.GetFaceVarArr(bound);
                var ft = fluxBound.FixT.GetFaceVarArr(bound);
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    double gamma = gf.GetFaceVar(i1);
                    double h = fh.GetFaceVar(i1);
                    double flux_a = fa.GetFaceVar(i1) - h;
                    double flux_b = fb.GetFaceVar(i1) + h * ft.GetFaceVar(i1);
                    flux_a *= face.Area;//实际上为a*Sb
                    flux_b *= face.Area;
                    double beforPhib = face.gDiff * gamma - flux_a;
                    double rightTerm = -gamma * face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, BouX[i1]) + flux_b;
                    double A = face.gDiff * gamma / beforPhib;
                    double B = rightTerm / beforPhib;
                    factors[i1 - start] = (A, B);//phi_b=A*phi_C+B
                    BouFaceShare[i1] = flux_b + flux_a * B;
                    ABouH[i1] = flux_a * A;
                }
            }
        }
        protected void formABOnFluxWallBound_Test(FvBound bound, FluxWallBound fluxBound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var gf = Gamma.GetFaceVarArr(bound);
            int start = bound.FaceIdx.Start;
            if (fluxBound.IsConstSource())
            {
                var factors = fluxBound.BouFactors1;
                var fb = fluxBound.FixFlux_b.GetFaceVarArr(bound);
                for (; i1 < endp1; i1++)
                {
                    (_, factors[i1 - start]) = formAb_Bound_FluxCombined(i1, gf.GetFaceVar(i1), 0, fb.GetFaceVar(i1));
                }
            }
            else
            {
                var factors = fluxBound.BouFactors2;
                var fb = fluxBound.FixFlux_b.GetFaceVarArr(bound);
                var fa = fluxBound.FixFlux_a.GetFaceVarArr(bound);
                var fh = fluxBound.Fixh.GetFaceVarArr(bound);
                var ft = fluxBound.FixT.GetFaceVarArr(bound);
                for (; i1 < endp1; i1++)
                {
                    double h = fh.GetFaceVar(i1);
                    factors[i1 - start] = formAb_Bound_FluxCombined(i1, gf.GetFaceVar(i1), fa.GetFaceVar(i1) - h, fb.GetFaceVar(i1) + h * ft.GetFaceVar(i1));
                }
            }
            (double A, double B) formAb_Bound_FluxCombined(int faceID, double gamma, double flux_a, double flux_b)
            {
                return formAb_Bound_FluxCombined_v1(faceID, gamma, flux_a, flux_b);
            }
            (double A, double B) formAb_Bound_FluxCombined_v1(int faceID, double gamma, double flux_a, double flux_b)
            {//利用几何关系进行了松弛
                ref var face = ref BouFaces[faceID];
                flux_a *= face.Area;//实际上为a*Sb
                flux_b *= face.Area;
                double svv = Vector3d.DotProduct(face.Normal, face.FaceNOrth) / InnCells[face.OSideCell].Volumn;
                double beforPhib = (svv + face.gT) * gamma - flux_a;
                double rightTerm = (-Vector3d.DotProduct(Grad[face.OSideCell], face.FaceNOrth) + svv * (BouX[faceID] - InnX[face.OSideCell])) * gamma + flux_b;
                double A = (svv + face.gT) * gamma / beforPhib;
                double B = rightTerm / beforPhib;
                BouFaceShare[faceID] = flux_b + flux_a * B;
                ABouH[faceID] = flux_a * A;//此边界应该是flux_a*边界值+flux_b，边界值是内部该网格的函数，因此可进一步线性化提高收敛性
                return (A, B);
            }
            (double A, double B) formAb_Bound_FluxCombined_v2(int faceID, double gamma, double flux_a, double flux_b)
            {//利用几何关系进行了松弛
                ref var face = ref BouFaces[faceID];
                flux_a *= face.Area;//实际上为a*Sb
                flux_b *= face.Area;
                double svv = Vector3d.DotProduct(face.Normal, face.FaceNOrth) / InnCells[face.OSideCell].Volumn;
                double beforPhib = (svv + face.gT) * gamma - flux_a;
                double rightTerm = (-Vector3d.DotProduct(Grad[face.OSideCell], face.FaceNOrth) + svv * BouX[faceID]) * gamma + flux_b;
                double A = face.gT * gamma / beforPhib;
                double B = rightTerm / beforPhib;
                //factors[faceID - start] = (A, B);  //phi_b=A*phi_C+B
                BouFaceShare[faceID] = flux_b + flux_a * B;
                ABouH[faceID] = flux_a * A;//此边界应该是flux_a*边界值+flux_b，边界值是内部该网格的函数，因此可进一步线性化提高收敛性
                return (A, B);
            }
            (double A, double B) formAb_Bound_FluxCombined_v3(int faceID, double gamma, double flux_a, double flux_b)
            {//利用上一步PhiC进行了松弛
                ref var face = ref BouFaces[faceID];
                flux_a *= face.Area;//实际上为a*Sb
                flux_b *= face.Area;
                double beforPhib = face.gDiff * gamma - flux_a;
                double rightTerm = gamma * (-Vector3d.DotProduct(Grad[face.OSideCell], face.FaceNOrth) + BouX[i1] * (face.gDiff - face.gT)) + flux_b;
                double A = face.gT * gamma / beforPhib;
                double B = rightTerm / beforPhib;
                //factors[faceID - start] = (A, B);  //phi_b=A*phi_C+B
                BouFaceShare[faceID] = flux_b + flux_a * B;
                ABouH[faceID] = flux_a * A;
                return (A, B);
            }
            (double A, double B) formAb_Bound_FluxCombined_v4(int faceID, double gamma, double flux_a, double flux_b)
            {//只是用了OverCor，也没有发散
                ref var face = ref BouFaces[faceID];
                flux_a *= face.Area;//实际上为a*Sb
                flux_b *= face.Area;
                double beforPhib = face.gDiff * gamma - flux_a;
                double rightTerm = -gamma * face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, BouX[i1]) + flux_b;
                double A = face.gDiff * gamma / beforPhib;
                double B = rightTerm / beforPhib;
                //factors[faceID - start] = (A, B);//phi_b=A*phi_C+B
                BouFaceShare[faceID] = flux_b + flux_a * B;
                ABouH[faceID] = flux_a * A;
                return (A, B);
            }
        }
        public override IConnectFactor GetConnectWallFactor(FvBound bound)
        {
            return new WallFactor(this, bound);
        }
        public class WallFactor : FactorCalcBase, IConnectFactor
        {
            public WallFactor(EquationWithDiffTerm equationWithDiffTerm, FvBound bound)
            {
                SetBasic(equationWithDiffTerm, bound);
                var matrixBound = equationWithDiffTerm.EquationBounds[bound.ID] as CoupledWallBound;
                this.gf = equationWithDiffTerm.Gamma.GetFaceVarArr(bound);
                this.fb = matrixBound.FixFlux_b.GetFaceVarArr(bound);
                this.fa = matrixBound.FixFlux_a.GetFaceVarArr(bound);
            }
            IFaceVarArr<double> gf;
            IFaceVarArr<double> fb;
            IFaceVarArr<double> fa;
            public (double a11, double a12, double b1, double a21, double a22, double b2) GetFactor(int faceID)
            {
                ref var face = ref BouCFaces[faceID];
                double gamma = gf.GetFaceVar(faceID);
                double gammadiff = -gamma * face.gDiff;
                double fcor = gamma * face.CalNonOrthogonalityCorrectionOfBound(Grad, InnX, BouX[faceID]);
                return (-gammadiff, gammadiff, fcor,
                 gammadiff, -gammadiff - fa.GetFaceVar(faceID) * face.Area, fb.GetFaceVar(faceID) * face.Area - fcor);
            }
            public void SetFactor(int faceId, double A11, double B1)
            {
                BouA[faceId] = -A11;
                BouFaceShare[faceId] = B1;
            }
        }
    }
}